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In this work, we present a general, efficient, and provably robust represen- 
tation for intrinsic triangulations. These triangulations have emerged as a 
powerful tool for robust geometry processing of surface meshes, taking a 
low-quality mesh and retriangulating it with high-quality intrinsic trian- 
gles. However, existing representations either support only edge flips, or 
do not offer a robust procedure to recover the common subdivision, that 
is, how the intrinsic triangulation sits along the original surface. To build 
a general-purpose robust structure, we extend the framework of normal 
coordinates, which have been deeply studied in topology, as well as the more 
recent idea of roundabouts from geometry processing, to support a variety of 
mesh processing operations like vertex insertions, edge splits, etc. The basic 
idea is to store an integer per mesh edge counting the number of times a 
curve crosses that edge. We show that this paradigm offers a highly effective 
representation for intrinsic triangulations with strong robustness guaran- 
tees. The resulting data structure is general and efficient, while offering a 
guarantee of always encoding a valid subdivision. Among other things, this 
allows us to generate a high-quality intrinsic Delaunay refinement of all 
manifold meshes in the challenging Thingil0k dataset for the first time. This 
enables a broad class of existing surface geometry algorithms to be applied 
out-of-the-box to low-quality triangulations. 


CCS Concepts: e Mathematics of computing — Mesh generation. 


Additional Key Words and Phrases: remeshing, intrinsic triangulation, De- 
launay triangulation, discrete differential geometry 


1 INTRODUCTION AND RELATED WORK 


Geometric data plays a growing role in applications from computa- 
tional fabrication to to autonomous driving to augmented reality, but 
data in these applications is increasingly difficult to deal with due 
to the poor quality of meshes generated by non-expert users, or by 
algorithms targeted at visualization rather than mesh processing— 
there have hence been significant recent efforts to make geometric 
algorithms more robust [Zhou et al. 2016; Hu et al. 2018; Sellan et al. 
2019; Sawhney and Crane 2020]. One basic tool is to remesh the 
input to obtain higher-quality elements, but for this approach to 
work on difficult, near-degerate inputs, remeshing algorithms must 
themselves be extremely robust. Moreover, traditional approaches 
to remeshing based on vertex positions in R” must negotiate the 
trade-off between mesh size, the quality of mesh elements, and 
geometric approximation of the input domain. 


Intrinsic Triangulations. A promising idea is to approach geome- 
try processing from the intrinsic point of view: rather than consid- 
ering the embedding of the geometry in space, one focuses only on 
point-to-point distances along the surface, as encoded by the edge 
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Fig. 1. We extend the machinery of normal coordinates beyond just edge 
flips, to enable a broader set of local mesh operations such as vertex insertion. 
By doing so, we dramatically improve the robustness of algorithms like 
Delaunay refinement in the intrinsic setting. 


lengths of a triangulation. This perspective is quite natural for prob- 
lems in geometry processing and scientific computing, since many 
objects in these domains are themselves intrinsic—for instance, the 
Laplace-Beltrami operator, which appears in numerous algorithms 
and fundamental partial differential equations (PDEs). In this paper 
we consider so-called intrinsic triangulations, whose edges no longer 
need to be straight line segements in Euclidean space, but can in- 
stead be any straight or geodesic path across the input polyhedron 
(see Figure 2). This construction still captures the input geometry 
exactly, but provides a dramatically larger space of possible triangu- 
lations, lending enormous flexibility to geometric algorithms. For 
instance, it de-couples the quality of elements used for simulation 
from the elements used to describe the geometry, side-stepping the 
trade-off encountered in traditional meshing. Moreover, this perpec- 
tive enables the input polyhedral surface to serve as a background 
domain—analogous to the Euclidean plane in traditional computa- 
tional geometry—allowing one to “port” trusted algorithms from 
the plane to curved surfaces. Most importantly, by encapsulating all 
this machinery in an interface that resembles an ordinary mesh, one 
can provide robustness as a subroutine: rather than make existing 
algorithms more robust one by one, we can transform the input 


Fig. 2. Edges of an intrinsic triangulation are allowed to be geodesic paths 
along a surface (left). The faces of such a triangulation can be laid out in 
the plane as ordinary triangles (right). 


into an intrinsic triangulation, execute an ordinary (“non-robust”) 
algorithm, and then read off the reults in a variety of ways. 

A remaining impediment to making the intrinsic approach truly 
reliable is to develop data structures for intrinsic triangulations that 
provide all the expected operations from standard mesh processing, 
while simultaneously providing strong guarantees of correctness. 
The basic challenge is encoding the correpsondence between the 
input mesh T° and an intrinsic triangulation T! sitting atop it, so 
that data on one triangulation can be transferred to the other. The 
first such data structure was the overlay mesh of Fisher et al. [2006]. 
The overlay uses a halfedge mesh decorated with special vertex and 
edge attributes to maintain the common subdivision of T° and T}, ie, 
the polygon mesh obtained by “slicing up” the underlying surface 
along the edges of both T° and T!. This approach guarantees cor- 
rect connectivity, but ordinarily-local operations such as edge flips 
become non-local and expensive to evaluate—moreover, edge flips 
are the only operation supported by this data structure. Sharp et al. 
[2019a] instead encode the correspondence implicitly by storing 
so-called signposts at vertices, which give the direction and length of 
each intrinsic edge. This approach is somewhat complementary to 
the overlay mesh: local mesh operations are now cheap to evaluate, 
but the encoding of connectivity now depends on floating-point 
values, and is hence not guaranteed to be correct. For instance, when 
tracing out intrinsic edges small floating point errors can cause one 
to “miss” the target vertex. A key development here, however, was 
extending intrinsic triangulations to operations beyond edge flips. 


Integer-Based Encoding. In this paper, we introduce an integer- 
based data structure for intrinsic triangulations that offers the best 
of both worlds: an implicit encoding of correspondence that sup- 
ports fast local operations, but which is also guaranteed to correctly 
describe connectivity. Like the signposts, our data structure also 
supports a wide variety of local mesh operations (Section 3). As 
demonstrated in Section 5, we get dramatically improved robustness 
for difficult tasks, e.g., we achieve a 100% success rate for extracting a 
high-quality Delaunay refinement of low-quality input data. In turn, 
any algorithm that relies on a high-quality triangulation (e.g., for 


Fig. 3. We build on the idea of normal coordinates count how many times a 
curve crosses each edge of a triangulation. 


solving PDEs) can immediately benefit from this improved robust- 
ness. For instance, in Section 4 we observe improved robustness for 
computing geodesic distance, local parameterization, constructing 
geodesics, and finding smooth vector fields. 

The basic starting point for our data structure is the concept of 
normal coordinates from geometric topology (not to be confused with 
geodesic normal coordinates from Riemannian geometry). However, 
we must augment this construction in several ways in order to make 
it suitable for geometry processing. As detailed in Section 2.3, the 
basic idea of normal coordinates is to simply count how many times 
each edge of a triangulation is crossed by some curve (Figure 3). 
Such coordinates were originally developed to study not curves, 
but rather embeddings of surfaces in 3-manifolds [Kneser 1929; 
Haken 1961; Hass and Trnkova 2020], and subsequently appear 
in several places in mathematics (e.g., for studying the mapping 
class group [Farb and Margalit 2011]), including significant work 
on algorithms [Bell 2015, 2018; Schaefer et al. 2008]. In theoretical 
computer science, normal coordinates are also viewed as a means of 
“compressing” curves, e.g., the total number of bits required to store 
a long winding curve can be exponentially smaller than storing 
explicit segments along the curve [Erickson and Nayyeri 2013]. 

One challenge with using normal coordinates for geometry pro- 
cessing is that existing literature rarely considers operations beyond 
edge flips: the little that does considers only closed loops e.g. [Schae- 
fer et al. 2002, Section 5.4], whereas curves that terminate at vertices 
are absolutely essential for encoding triangulations. A second issue 
is that normal coordinates alone are not enough to uniquely identify 
curves implied by the coordinates with logical edges of a mesh. Very 
recently, Gillespie et al. [2021, Section 5.2] proposed a solution to 
this issue using what they call roundabouts, but again do not con- 
sider operations beyond edge flips. Third, whereas most literature 
assumes that normal coordinates encode homotopy classes of curves 
in a purely topological setting (or perhaps hyperbolic geodesics), 
we must make a significant departure this perspective and assume 
that the normal coordinates encode a triangulation of a Euclidean 
polyhedron by geodesic edges. This distinction is important since, 
in general, not all normal coordinates describe a valid Euclidean 
geodesic triangulation. Considering this special case in turn enables 
us to establish procedures not previously seen. 


Contributions. Overall, we make the following contributions: 


e We describe normal coordinates as a representation for gen- 
eral intrinsic triangulations, including the case where vertices 
have been added to the triangulation. 

e We extend integer-based data structures for geometric intrin- 
sic triangulations to include local operations beyond edge 
flips. 

e We prove the correctness and quality of a Delaunay refine- 
ment algorithm for intrinsic triangulations of surfaces with- 
out boundary. 

e We also extend intrinsic Delaunay refinement to surfaces 
with boundary. 

e We introduce a new, more accurate way of transferring func- 
tions between bases on different triangulations. 


We also experimentally validate the robustness of our technique, 
including generating intrinsic Delaunay refinements for all manifold 
meshes in the Thingi10k dataset, and demonstrate robustness for a 
variety of basic algorithms from geometry processing. 


2 NOTATION AND CONVENTIONS 
2.1 Connectivity 


Throughout we assume that our domain is an oriented 
manifold surface M, possibly with boundary. We write 
T = (V,E, F) to denote a triangulation of M with ver- 
tices V, edges E, and faces F. In general we allow tri- j >j 
angulations that are not be simplicial, but can instead 

be a A-complex in the sense of Hatcher [2002, Section 

2.1]—we allow, e.g., two edges of the same triangle to 

be glued together (see inset). We will refer to vertices jj 

i € V, edges ij € E, and faces ijk € F by one, two, or three indices, 
resp.. Note that as our triangulations need not be simplicial, the 
vertices i, j of an edge may not be distinct, and do not necessarily 
identify the edge—there may be multiple edges between i and j. 
We will refer to oriented halfedges ij € H, and will 
use ul K to denote a value u at corner i of f triangle ijk, 
and ujj to denote a value u at halfedge ij. Through- 
out, we consider a fixed input triangulation T° of 
M, as well as a dynamic intrinsic triangulation T! sitting atop M. 
T! must contain all vertices of T°, ie. V! D V?, but may include 
additional inserted vertices which we denote by V* := V! \ V°. We 
will generally use indices a, b, c for vertices of T° and indices i, j, k 
for vertices of T! (which may also be in T°). 


2.2 Geometry 


The geometry of a triangulation is determined by a collection of 
edge lengths £ : E — Rs» satisfying the triangle inequalities in each 
face. For instance, we typically begin by assigning E° and E! the 
same edge lengths €; ; £i; = |fj — fil determined by vertex coordinates 
f: V —> RÌ. While ‘the lengths £? remain fixed, the edge lengths 
t! of T! may change due to operations like intrinsic edge flips 
(Section 3.3). Even after changing the edge lengths, each individual 
triangle ijk € F! can always be drawn as an ordinary triangle in 
the Euclidean plane, allowing us to compute quantities such as face 
areas or corner angles. In general, we will not need to simultaneously 
embed all triangles of T! in RÌ. Finally, we use exp, : TyM — M 
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Fig. 4. Points are encoded relative to both triangulations T° and T!. For each 
triangulation we store the simplex containing the point, and the barycentric 
coordinates within that simplex. Here we show a few examples. 


to denote the exponential map at x. Given a 
tangent vector u at x, expx(u) is the point 
reached by walking straight along the surface 
in the direction of u for a distance |u| (see 
inset). (In practice this can be implemented as 
in [Sharp et al. 2019a, Section 3.2.2]). 

Points x € M can be expressed in barycen- 
tric coordinates relative to some simplex (ver- 
tex, edge, or triangle) of a triangulation T, 
e.g., a point in a triangle ijk is given by 
three coordinates uj, uj, uk € [0,1] such that 
uj + uj + ug = 1. For vertices, we set the sin- 
gle coordinate u; to 1. Importantly, we will 
need to encode points x with respect to two 
different triangulations T° and T!, using barycentric coordinates 
u and v, resp. (Figure 4). Note that for vertices that are shared by 
both triangulations, we have uj = vj = 1. We will use q? to denote 
a location on T° represented by simplex along with a barycentric 
coordinate, and similarly will use q! for a location on T!. 


exp,(u) 
inas 


2.3 Integer Coordinates 


Normal coordinates. We use normal coordinates to count the num- 
ber times each edge of T! crosses the edges of T° (see Figure 5, left). 
In principle, normal coordinates could either be defined as a value 
per edge of T°, counting the number of crossings from T!, or as a 
value per edge of T!. In our setting, we take the latter approach, as it 
remains fully-informative even when we insert new vertices in to T1. 
It is then natural to think of T! as the primary triangulation, with 
T° as a collection of geodesic curves sitting along it (Figure 5, right). 
To emphasize this abstract viewpoint, we will refer to the edges 
of E° as curves. Because we allow edges to be split, a single edge 
ab € E? may actually correspond to a sequence of curves expressed 
in normal coordinates which meet at intermediate vertices i € V*. 

Precisely, we use n : E! — Z to denote the normal coordinates. 
For each edge ij € E!, the quantity nj; indicates how many times 
ij is crossed by edges of E°: if nij > 0, then this is a count of 
crossings (formally, transversal intersections), whereas if nj; = —1, 
it indicates that a curve runs along edge ij (Figure 5, left). Note that 
nij is never less than —1 because we are working with triangulations, 
and multiple edges of a triangulation cannot lie along the exact same 
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Fig. 5. The intrinsic triangulation T! is often presented embedded in R? on 
top of the input triangulation T° (left). However, in our setting it is helpful to 
think of T! as an abstract intrinsic triangulation, given only by connectivity 
and edge lengths, which carries a collection of geodesic curves (right). 


path. We use nj, to denote the number of transversal crossings, i.e. 
ni; := max(nij, 0), and similarly define ni; := —min(njj,0), which 
is 1 on shared edges and 0 otherwise. 

Additionally, we define two quantities at each corner of T!, which 
count how many curves emanate from the interior of the corner 


and cross the corner resp. (Figure 6): 


p 

e) := max(0, nij — ny - ngi) (1) 
ij 1 + + _ +) _ ik _ ki 

C= 3 (max (0. nat Nki nt) ee; ) : (2) 


Crossings. A point where a curve crosses an edge ij € E! can 
be described either in a combinatorial or geometric sense. A com- 
binatorial crossing is given by a pair ¢ = (ij, p), such that ¢ is the 
ph crossing along oriented edge ij € H!. A geometric crossing is 
similarly given by z = (ij, p, u,v), where u, v encode the location of 
the crossing along the curve and along the ij (resp.) in barycentric 
coordinates. Importantly, both of these crossings are oriented: the 
choice of halfedge ij versus ji indicates which side of the edge 
one is “coming from” and “going to? which will be important when 
tracing out curves along the surface. We let ¢ := (ji, nij -p-1) 
denote a reversal of orientation. 


Roundabouts. Normal coordinates alone do not fully encode the 
correspondence between T° and T! because we cannot necessarily 
determine which curve along T! corresponds to which edge of T° 
(recall that there may be multiple edges of T° between the same pair 
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edges leaving corner k edges crossing corner k 


Fig. 6. It is often useful to count the curves emanating from the interior of 
each corner (left), and the curves crossing each corner (right). 
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Fig. 7. We employ roundabouts to encode how edges of T° and T! are 
interleaved around vertices. 


of vertices). Thus we additionally store roundabouts r : H! > Z>0; 
introduced by Gillespie et al. [2021, Section 5.2], which describe 
how the edges of T° and T! are interleaved around vertices. Unlike 
Gillespie et al., we may insert new vertices—however, we still store 
roundabouts only at halfedges pointing away from shared vertices 
a € V’. These are sufficient to disambiguate the identity of any 
traced curve, since all edges of E° must start and end at vertices in 
V° (if an edge has been split, then the sequence of curves starts and 
ends at vertices in V°). 

Precisely, for each halfedge aj € H! starting at a shared vertex 
a € V°, the roundabout stores the first halfedge ab € H? following 
aj. This is encoded as an index rgj € Z>0, where we enumerate the 
halfedges of T° about vertex a in counterclockwise order (Figure 7) 


3 ALGORITHMS AND DATA STRUCTURES 


In this section, we provide descriptions of our data structure and 
the operations that it supports. Detailed pseudocode can be found 
in Appendix A. 


3.1 Data Structure 


The most essential data to maintain is a mesh of triangulation T! = 
(V!,E!, F!); recall that V! > V°. We use a halfedge mesh, since 
halfedge meshes can represent general A-complexes (Section 2.1), 
though one could also use a vertex-face adjacency list plus a small 
amount of additional data [Sharp and Crane 2020a, Section 4.1]. On 
top of this mesh, our data structure maintains four quantities: 


e lengths £i; € Ryo for each edge ij € Et, 

e normal coordinates nj; € Z for each edge ij € E}, 

e roundabouts raj € Z>o for each halfedge aj € H! incident 
on a vertex a shared with T°, 

e barycentric coordinates q? relative to T° for each i € V1. 


This differs from the scheme of Gillespie et al. [2021, Section 
5] in a few key ways. The essential difference is that their data 
structure assumes that T° and T! share the same vertex set (V! = 
V°). This has numerous consequences—e.g. they use nonnegative 
normal coordinates nij € Z>o0, they do not store input positions 
q. Most importantly, it is impossible to perform many of the local 
mesh operations we describe in the next section without changing 
the vertex set of T1; our generalization of the representation is 
essential if one wants to perform tasks such as Delaunay refinement 
(Section 4.2). 


Case 2 


Fig. 8. A curve entering triangle jik along edge ij can proceed in 3 ways: it 
can exit along edge ik, in which case it is counted by ck (left); it can exit 
along edge kj, in which case it is counted by cik (center); or it can terminate 
at vertex k (right). This forms the core of procedure TRACEFROM. 


3.2 Extracting Curves 


Our first task is to recover a curve on T! from its normal coordinates— 
because our curves are geodesic, we can determine the exact ge- 
ometry from these normal coordinates. In particular, we describe a 
procedure EXTRACTCURVE (Algorithm 2) which takes in any combi- 
natorial crossing ¢ along a curve, and computes the curve’s trajec- 
tory along T! as a sequence (i, z1,..., Zk, j) of geometric crossings 
along with start and end vertices i, j € V1. We note that this mirrors 
the discussion in Gillespie et al. [2021, Section 6], albeit in a more 
general setting; we include a full description here for completeness. 

We proceed in two steps, first determining the triangle strip that 
the curve passes through and only then computing the curve’s 
geometry. Note that the triangle strip depends solely on the integer- 
valued normal coordinates n, while geometric data and floating 
point computation are relegated to the second step. 

The first step is performed by TRAcEFRom (Algorithm 1), which 
takes some combinatorial crossing ¢ = (ij, p) along a curve y, and 
traces out the remaining combinatorial crossings until y terminates 
at a vertex. TRACEFROM proceeds iteratively, taking the crossing 
where y enters a triangle, and using the triangle’s normal coordinates 
to determine where y exits (see Figure 8). The direction in which to 
trace the curve is determined by the orientation of ij. 

To determine the triangle strip containing y, ExTRACTCURVE calls 
TRACEFROM once in either direction, yielding the sequence of all 
combinatorial crossings along y. ExTRacTCuRVE then unfolds this 
triangle strip in an arbitrary planar coordinate system and draws 
y as a straight line between its endpoints. The intersection of this 
line with each of the intermediate edges determines the geometric 
crossings along y (Figure 9). Note that unlike Gillespie et al. [2021, 
Algorithm 1], we may have to invoke ExrracTCurVE multiple times 
to extract a single edge ab € E°, as it may pass through several 
vertices of T1, e.g. due to edge splits (Section 3.5). 

Finally, it is sometimes useful to convert a single combinatorial 
crossing ¢ into a geometric crossing z. We will refer to this operation 
as EXTRACTGEOMETRICCROSSING; it may be implemented by calling 
EXTRACTCuRVE and then returning the single desired crossing. 


3.3 Edge Flip 


Edge flips are a well-studied operation, but we include a discussion 
here for completeness. We may flip and edge if and only if (i) both 
endpoints have degree at least one after the flip, and (ii) the two 
triangles containing the edge form a convex quadrilateral (Figure 10). 


Fig. 9. We compute barycentric coordinates by laying out a triangle strip in 
the plane. 
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Fig. 10. An edge flip replaces an edge with its opposite diagonal (left). An 
edge ij is not be flippable if it would leave a vertex with degree zero, or if 
its neighboring faces form a nonconvex quadrilateral (right). 


Mesh Update. We replace edge ij with an edge kl. We compute 
the new edge length th by laying out the two old triangles ijk, lji 
in the plane and measuring the length of the appropriate diagonal. 


Normal Coordinates & Roundabouts. The new normal coordinate 
nj does not depend at all on the geometry of T!. It is given by the 
following formula: 


— Jk, t [il _ oki 
nki = Cj +c ts Cj cj 


el tei telte an 
i i J J 


1 


Lj 
e — © 


Ci L 


jk|_ 1J _ 1 
i | Ie 2. 
This differs slightly from the formula of Gillespie et al. [2021], which 
did not allow for inserted vertices. 

We can update each roundabout from its previous neighbor: 


rei = mod (rg + ej! + nz, dego(k)), (4) 


where deg,(k) is the degree of vertex k in triangulation T°. The 
quantity eil + nz; counts how many edges of T? are between ki and 
kl: eil counts edges strictly between 1 them, and nọ; adds one if there 
is also an edge lying exactly along kl. 


We only perform this update for halfedges whose source is in V°. 


3.4 Face Split 


We now describe procedure SpLiTFAce (Algorithms 4 and 5), the first 
of several new routines to mutate triangulation T! while tracking the 
correspondence with T°. Note that Schaefer et al. [2002, Section 5.4] 
describe a similar face split operation in the topological setting, but 
do not provide the ability to insert a point at a particular geometric 
location, which is essential in our setting of Euclidean polyhedra. 
In particular, suppose we wish to insert a vertex at a point x € M, 
given by barycentric coordinates v on a triangle ijk € F!. To do 
so, we need to update the connectivity T1, edge lengths ¢!, normal 
coordinates n, and roundabouts r. Additionally, we need to compute 
the position q? of this new vertex in barycentric coordinates on T°. 


Mesh Update. We update the connectivity of T! with a new vertex 
and three new edges and faces. We compute new edge lengths as a 
formula of the barycentric coordinates u. Schindler and Chen [2012, 
Section 3.2] show that the length of a displacement vector du; in 
barycentric coordinates is give by 


\|Suil|? = -€7,duiSuj — pôu juz — h Ôupôui. (5) 


Normal Coordinates & Roundabouts. Unlike the 
case of an edge flip, where the new normal coor- 
dinate depends solely on the initial normal coordi- 
nates, vertex insertion is an inherently geometric 
operation. Different points necessarily result in dif- 
ferent normal coordinates (see inset), depending on 
the region R in which the point lies. 

Concretely, we first compute the geometric cross- 
; ings of all curves passing through face ijk (using 
the EXTRACTGEOMETRICCROSSING subroutine). We 
then determine which region R the new point lies 
in via a series of line-side tests. (One might in prin- 
ciple be able to reduce the number of curves that 
need to be extracted via lazy evaluation or caching, 
though we do not pursue such optimizations here.) 
We may misclassify points extremely close to a 
, region’s boundary due to floating point error, in 

which case we insert a valid point in the identified 
region, at a virtually identical location. Note that this behavior is per- 
fectly reasonable in, e.g., retriangulation algorithms (see Section 4.2), 
where the insertion location is not computed exactly anyway. 

The roundabouts on any new halfedges emanating from original 
vertices (i.e., vertices in {i, j, k}n V?) can be set from their neighbors 


via Equation 4. 


Position on T°. To determine q°, we must locate the triangle 
abc € F? containing x, as well as the barycentric coordinates of x 
within abc. We do so via interpolation from the corners of R. Ex- 
plicitly, the corners of R are all geometric crossings with known 
barycentric coordinates in some triangle abe € F° (computed in 
EXTRACTGEOMETRICCROSSING); we can then solve a small linear 
system to recover the barycentric coordinates of p in the same trian- 
gle. Intuitively, we recover generalized barycentric coordinates for 
p with respect to the polygon R and apply them on abc to recover 
standard barycentric coordinates in T° (see Appendix B for details). 


SpLiTFace(ijk, (uj, uj, 0)) 
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Fig. 11. Generally, one can split edge ij by performing a face split on a 
neighboring face followed by an edge flip (top). However, if ij carries a 
curve, this strategy will cause the insertec vertex to miss the curve (bottom). 
We hence provide a different edge split procedure for this case in Section 3.5. 


3.5 Edge Split 

We also introduce an operation SPLITEDGE (Algorithm 6), which 
takes as input a point given by barycentric coordinates v along an 
oriented edge ij € H!. If ij does not have a curve running along it 
(i.e. nij = 0), then this is implemented as a face split followed by an 
edge flip (Figure 11, top). However, if njj < 0 (which is common in 
practice—e.g. Section 4.2), we perform an explicit edge split which 
inserts the new vertex along the coincident curve (Figure 11, bottom). 


Mesh Update. We insert a new vertex and triangulate any adjacent 
faces, computing the new edge lengths via Equation 5. 


j Normal Coordinates & Roundabouts. When nij <0 
the new normal coordinates are simple functions of the 
old ones, since every curve in face ijk must emanate 
P from i or j, or cross k. The number of such curves 
is max(ng;, jx, 0) and edge pk crosses them all. We 
hence set np; and np; equal to nj;, and set 


i npk = Max(Np;,N jk» 0). (6) 


As with face splits, roundabouts on any new halfedges emanating 
from original vertices can be set from their neighbors (Equation 4). 


Position on T°. To determine q°, we must locate the edge ab € E? 
containing x, as well as the barycentric coordinates of x within ab. 
Since, i and j necessarily have known locations along some edge in 
ab € E°, we can simply interpolate by v to compute the location q°. 


3.6 Vertex Removal 


In general, a vertex which is present in the original triangulation 
cannot be removed without distorting the intrinsic metric because 
any curvature at that vertex would be lost. However, inserted ver- 
tices i € V* have no curvature, and can hence be removed safely. 
In fact this operation will be necessary for Delaunay refinement of 
domains with boundary (Section 4.2). 


2 33 Case 2 


Fig. 12. We extract the connectivity of common subdivision within each 
triangle using its normal coordinates. 


The basic strategy behind REMOovEVERTEX (Algorithm 7) is to 
flip edges incident on the vertex to be removed until it has degree 
three, then delete the three edges incident on the vertex as well 
as the vertex itself. No other data needs to be updated, since the 
edges of the resulting triangle already appear in the triangulation. 
Algorithm 7 describes this procedure, and Theorem D.1 proves its 
correctness for simplicial complexes. A nearly identical procedure 
can be used to remove an inserted boundary vertex. Schaefer et al. 
[2002, Section 5.4] also suggest a similar flipping procedure, but work 
in the topological setting where the necessary edge flips are always 
valid—they do not consider the convexity condition (Section 3.3). 


3.7 Moving Inserted Vertices 


Given the previous operations, we can easily define a procedure for 
moving around inserted vertices. Specifically, given a vector v in 
the tangent space of an inserted vertex i, we can move i along v in 
the following way: 


e First, compute the new location p = exp;(v). 
e Insert p using SPLITFACE. 
e Remove i using REMOVEVERTEX. 


We insert p first since the removal procedure could flip edges in- 
cident on the triangle containing p, invalidating its barycentric 
coordinates. Note that Sharp et al. propose an alternative strategy 
for local vertex displacement [Sharp et al. 2019a, Section 3.3.3]. 


3.8 Common Subdivision 


As noted previously, the common subdivision S of T? and T! is the 
polygon mesh obtained by “slicing up” the underlying surface along 
the edges of both T° and T!. The vertices of S are hence a superset 
of V? and V!, and every edge or face of T? and T! can be expressed 
as union of edges or faces of S (resp.). Moreover, the faces of S are 
always planar and convex. Most importantly in our setting, any 
piecewise-linear function on T° or T! can be represented exactly as 
a piecewise-linear function on S. Note however that even if T° and 
T! have nice elements, S is not in general a high-quality mesh, and 
may not itself be suitable for, e.g., solving PDEs. Rather, it plays a 
complementary role in the geometry processing pipeline, enabling 
(for instance) transfer of data between triangulations (Section 4.4), 
or visualization of data downstream via standard rendering tools. 


We compute the common subdivision by cutting T! along the 
edges of T°. First we extract the connectivity of S, using only the 
normal coordinates n;;. Then we recover the intersection geometry, 
allowing us to interpolate data stored at the vertices of T? or T! to 
S—most commonly, vertex positions on T° along with any solution 
data on T!. Note that this procedure was previously described by 
Sharp et al. [2019a, Section 3.4.2]; we recap it here for completeness, 
and to give a convenient description using our integer coordinates. 


Connectivity. We subdivide T! independently in each face ijk. The 
normal coordinates alone determine the connectivity of S within this 
face. There are just two cases to consider, illustrated in Figure 12. 
Case 1 occurs when no curves emanate from any corner, so we 
simply need to connect the first clk crossings along edge ij to the 
first a” crossings along ik (in order), and likewise for corners j 
and k. In Case 2 curves emanate from some corner; without loss of 
generality, let this corner be k so that the number of such curves is 
e! > 0 and there are more curves crossing edge ij than the other 
two edges. Hence, we can walk from i to j, connecting the first c k 
crossings to those along ik, the next elk crossings to vertex k, and 
the remaining ki crossings to those along edge k j. Note that curves 
running along edges (nij < 0) require no special treatment. 


Intersection Geometry. Next, we associate each vertex i of the 
common subdivision with a point in T° and a point in Tt, encoded 
in barycentric coordinates relative to some simplex (vertex, edge, or 
face) of the appropriate triangulation. Using this, one can linearly 
interpolate data in the usual way. Again, there are just two cases: 
each vertex i in S is either a vertex of T! or the intersection of an 
edge of T° with an edge of T!. In the first case, the position on T! is 
given by i itself, and its position q? on T° was computed when i was 
inserted. In the second case, we compute the desired barycentric 
coordinates using ExTRACTCURVE (see Algorithm 8 for details). 


3.9 Transposing Coordinates 


Throughout, we store normal coordinates n : E! > Z which count 
how many times edges of E! cross edges of E°. It is sometimes 
useful to observe that tracing an edge ab € E? over T! counts how 
many times ab crosses edges of T!. If T! has more vertices than T°, 
then the edges of E! are not normal over T°—they can start and 
end in the middle of faces of T°—and these crossing counts do not 
uniquely encode the structure of T1. However, if T° and T! do have 
the same vertex set, then these crossing counts provide an implicit 
representation of T! as a collection of curves over T°. 


3.10 Visualization 


In the following, we show examples of intrinsic triangulations on top 
of meshes (e.g. Figure 13). To produce these figures, we compute the 
common subdivision (Section 3.8) and draw the edges of the input 
mesh with a black wireframe while coloring the intrinsic triangles in 
arbitrarily-chosen colors. In figures displaying functions defined on 
intrinsic triangulations (e.g. Figure 15), we interpolate the solutions 
along the common subdivision for rendering. 


3.11 Robust Implementation 


Our integer coordinates are guaranteed to encode a triangulation 
sitting atop T!. The geometric accuracy of this triangulation, of 
course, depends on floating point arithmetic, which can become 
inaccurate in near-degenerate configurations. Exact predicates have 
been applied with great success to similar problems [Devillers and 
Pion 2003]. Unfortunately they do not directly apply to intrinsic 
triangulations, as the predicates that we evaluate are not fixed func- 
tions of the input data; an intrinsic edge length can depend upon 
arbitrarily many input edge lengths. Hence, we focus on fast and 
robust implementations using ordinary floating point arithmetic. 
One essential tool for dealing with intrinsic triangulations on 
near-degenerate input meshes is intrinsic mollification, introduced 
by Sharp and Crane [2020a]. Mollification improves degenerate 
meshes by adding a small € to every edge length, provably improving 
triangle quality. This changes the geometry by a negligible amount, 
and moreover we only mollify if some triangle is within e of being 
degenerate. This procedure works particularly well with our data 
structure compared to signposts: the signpost data structure relies 
on tracing queries along the surface which become less accurate 
when mollification is applied. Our integer coordinates have no such 
problem: we always get the correct edge sequence, even if the mesh 
geometry is slightly modified. In our experiments we mollify with 
e = 107°, and find that it resolves almost all numerical difficulties. 
Even after mollification, it is still beneficial to use care when 
working with floating point. For example, there are well-conditioned 
triangles on which the Delaunay condition (Equation 8, discussed in 
the next section) is difficult to evaluate; in practice, we only enforce 
Equation 8 up to some e tolerance. As a further example, when 
computing new normal coordinates in SpLiTFace, one could lay 
out the face in the plane, and independently count intersections 
along the new edges. However, this can produce invalid normal 
coordinates in floating point. We apply a more complicated policy 
(see Appendix A) which always yields valid normal coordinates. 
For additional details on all procedures, we refer the reader to our 
implementation, which will be made available after review. 


3.12 Other Algorithms 


Normal coordinates also enable a wide variety of other operations 
not detailed here. For instance, Schaefer et al. [2002, Section 5] pro- 
vide algorithms for counting connected components, checking if 
crossings are part of the same curve, checking if curves are iso- 
topic, and computing the oriented intersection number. Erickson 
and Nayyeri [2013] provide an asymptotically-fast algorithm for 
tracing normal curves across a surface. Finally, Dynnikov [2020, 
Proposition 13] provides an algorithm for computing how many 
times curves represented by normal coordinates intersect. 


4 APPLICATIONS 
4.1 Intrinsic Delaunay Triangulations 


One key application of intrinsic triangulations 
is the computation of intrinsic Delaunay trian- ‘ 
gulations (Figure 13, top). A triangulation is said 
to be Delaunay if the sum of angles opposite 


Fig. 13. Using our integer-based data structure, we can not only improve 
near-degenerate meshes by generating intrinsic Delaunay triangulations 
(top), but can also extract the common subdivision after computing a high- 
quality intrinsic Delaunay refinement (bottom). 


every edge is at most 7, i.e. for every ij € E we have 
a! + oy <r. (7) 


Delaunay triangulation have a number of beneficial properties. 
One consequence of Equation 7 is that edges of a Delaunay triangu- 
lation must have nonnegative cotan weights: 


cot 0 + cot 0f > 0. (8) 


In fact, Equation 8 is equivalent to Equation 7 above, and provides 
a convenient formula for checking the Delaunay property in an 
intrinsic triangulation. Moreover, Equation 8 ensures that the finite 
element Laplacian L satisfies the maximum principle, guaranteeing 
that discrete harmonic functions do not have local extrema in the 
interior of the domain [Bobenko and Springborn 2007, Proposition 
19]. Similarly Equation 8 also ensures that discrete harmonic vector 
fields are “flip-free” [Sharp et al. 2019b, Section 5.4]. Furthermore, 
the local Delaunay condition implies the empty circumcircle property: 
each triangle’s geodesic circumdisk contains no vertices, illustrated 
in Figure 14, left [Bobenko and Springborn 2007, Proposition 10]. 
The Delaunay triangulation can be computed via a simple greedy 
algorithm: flip any non-Delaunay edge until all edges satisfy Equa- 
tion 7 [Bobenko and Springborn 2007, Propositions 11 and 12]. 


Fig. 14. Triangles in Delaunay meshes have empty circumdisks, and thus 
well-defined circumcenters (left). When necessary, we locate a triangle’s 
circumcenter by walking outwards from its barycenter (right). 


4.2 Intrinsic Delaunay Refinement 


Delaunay refinement inserts vertices in order to produce a Delaunay 
mesh whose triangles all satisfy a minimum angle bound (Figure 13, 
bottom). Here we modify Chew’s second algorithm to perform in- 
trinsic Delaunay refinement [Chew 1993; Shewchuk 1997]. This 
problem has been extensively studied in the plane, but an intrinsic 
(i.e. geodesic) scheme was only recently proposed by Sharp et al. 
[2019a, Section 4.2]. However, they did not handle meshes with 
boundary—here we resolve the essential difficulties of the bound- 
ary case, and show how refinement can be implemented using our 
integer-based data structure. 

In the plane, the basic algorithm is to greedily pick any triangle 
which violates the minimum angle bound, insert a vertex at its cir- 
cumcenter, then flip to Delaunay. This process continues until all tri- 
angles satisfy the angle bound. If a triangle’s circumcenter is outside 
the domain, then the boundary edge ij separating the triangle from 
its circumcenter is split at its midpoint; subsequently, all interior 
vertices within at least a distance of £;;/2 are removed—though re- 
moving additional interior vertices causes no issues (Appendix C.1). 
One can prove that this process succeeds for minimum angle bounds 
up to 25.65 degrees on planar domains with boundary angles at least 
60° [Shewchuk 1997, Section 3.4.2]. More advanced versions of this 
procedure can achieve better angle bounds, e.g. [Rand 2011], but 
here we restrict our attention to the basic algorithm for simplicity. 

There are two difficulties in adapting this algorithm to the in- 
trinsic setting: locating circumcenters and computing (geodesic) 
distances. As mentioned earlier, intrinsic Delaunay triangulations 
obey the empty circumcircle property; hence each triangle has an 
intrinsically-flat circumdisk with a well-defined center (Figure 14, 
left). So long as this center corresponds to a point on the surface, it 
can be found by walking from the triangle’s barycenter (Figure 14, 
right). In practice, we compute triangle ijk’s circumcenter in ho- 
mogeneous (i.e., unnormalized) barycentric coordinates ô; via the 
following formula [Schindler and Chen 2012, Section 2.3]: 


> sigh gh gh ap 
ôi = CE, + ey — Gs (9) 
and then normalize to obtain barycentric coordinates 
.— Îi 
Uj := O40) 40K” (10) 


To locate the circumcenter on the surface, we then evaluate the 
exponential map (Section 2.2) starting at the barycenter w; = wj = 
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Fig. 15. Running PDE-based algorithms such as the heat method on poor 
triangulations (left) can lead to inaccurate solutions. Flipping to intrinsic De- 
launay (center) and performing Delaunay refinement (right) can drastically 
improve the results. 


wk = 1/3, along the vector v — w. If we hit a boundary edge ij while 
tracing out this path, then the circumcenter is not contained in the 
surface, so we split ij at its midpoint and flip to Delaunay. We must 
then remove all inserted interior vertices within a geodesic ball 
of radius ¢;;/2 centered at the inserted point. Computing geodesic 
distance on a surface mesh is nontrivial, but Xia [2013, Corollary 1] 
shows that on a Delaunay triangulation any vertex inside a geodesic 
ball of radius r will also be inside the Dijkstra ball of radius 2r (i.e. 
points whose distance along the edge graph are at most 2r). We 
hence remove all interior inserted vertices within a Dijkstra distance 
of fij. Note that while Xia considers only the planar setting, their 
proof (which is based on triangle strips) applies without modification 
to intrinsic Delaunay triangulations of surfaces. 
On meshes with narrow cone vertices or boundary 60° 

angles, it may be impossible to find any triangulation 

satisfying a given angle bound. In such cases, we do 

not insert circumcenters of intrinsic triangles which are 

incident on exactly one narrow vertex, or are entirely 

contained in a triangle of T° which is incident a narrow 

vertex, and ignore such triangles when computing the 

minimum corner angle of the output mesh. Although 

the final output may violate the angle bound, such trian- 

gles appear only near narrow vertices. In analogy with 

the planar case, we set 60° as the minimum allowed 

angle sum (see inset); in practice the vast majority of 

meshes obey this constraint at all vertices (97.2% of Thingil0k), and 
even on those which do not we obtain high-quality triangulations. 


4.3 PDE-Based Geometry Processing 


PDE-based methods abound in geometry processing, as they are 
generally simple to implement and benefit from decades of research 
into linear solvers, Many such methods depend only on intrinsic 
data, and are hence a natural application of our intrinsic Delaunay 
triangulations and refinements. On near-degenerate inputs, simply 
running the standard algorithm on an intrinsic triangulation instead 
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Fig. 16. Here we compute a local parameterization (the logarithmic map, top), 
and a smooth vector field (bottom) using the connection Laplacian. Both pro- 
cedures yield inaccurate results on near-degenerate inputs (left)—intrinsic 
Delaunay triangulations (center) and intrinsic Delaunay refinements (right) 
greatly improve solution quality. Whether our solution is a scalar function 
or vector field, we can visualize it on the common refinement. 


of the original mesh yields solutions of dramatically higher quality 
( Figures 15 and 16). 

We show several examples: fast geodesic distance computation 
[Crane et al. 2017], local parameterization via the logarithmic map 
[Sharp et al. 2019b], and smooth vector fields [Knöppel et al. 2013]. 
Further examples on tasks such as parameterization, minimal sur- 
faces, and surface editing can be found in [Sharp et al. 2019a; Sharp 
and Crane 2020a]. Across the board, normal coordinates offer im- 
proved robustness guarantees, and open doors to higher solution 
accuracy with the common subdivision. 


4.4 Attribute Transfer 


Intrinsic triangulations can drastically improve the quality of so- 
lutions to PDEs on low-quality meshes. However in practice, one 
often needs to represent the solution on the input mesh. Past ap- 
proaches have simply “copied back” the solution values at vertices 
of the original mesh, but this strategy is ad-hoc and suboptimal. A 
more principled approach is to choose the function on the original 
mesh which is closest to the intrinsic solution. Here, we restrict 
our treatment to piecewise-linear bases and L? distance for sim- 
plicity, though the same strategy could easily be applied to other 
basis functions and notions of distance. Precisely, given a function 
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Fig. 17. Accuracy is improved by transferring PDE solutions back to an 
original triangulation as the L?-nearest solution, evaluated via the com- 
mon subdivision. Here we generate random low-quality meshes of the unit 
square by random edge splits (left), and plot the error in the solution of a 
Poisson equation compared to analytic ground truth, always represented in 
the basis of the original triangulation (right). Each data point is the average 
error over 100 trials. As expected, solving on the intrinsic Delaunay triangu- 
lation dramatically increases accuracy, but further improvements are gained 
by choosing the solution on the original mesh which is L?-nearest to the 
intrinsic solution, rather than naively copying vertex values. 


f on the intrinsic triangulation, we seek Ô on the original mesh that 
minimizes the squared L? distance 


If- ÎI? = J KOROKA (11) 


Here, f and Ô are functions represented in finite-dimensional bases 
with nodal values at the vertices of the intrinsic triangulation and 
the original mesh resp. In traditional finite elements, this integral 
commonly arises over a single triangulation, in which case it can be 
evaluated via the Galerkin mass matrx M as 


== =M - f), (12) 
where M is constructed as in [Strang and Fix 2008, Chapter 10, 
(32)]. However, in our setting f and f are encoded over different 
triangulations; they are members of different function spaces. Our 
key observation is that the common subdivision S (Section 3.8) 


provides exactly the structure needed to evaluate || f — || on as both 
functions are linear on each triangle of S. In fact, we have 


If- ÎI} = (Pf - Pof)” Ms (Pif - Pof), (13) 


where now Ms is the Galerkin mass matrix of the common subdi- 
vision, and Po, P; are interpolation matrices which map piecewise- 
linear functions on original and intrinsic triangulations to piecewise- 
linear functions on S, resp. In particular, Po is a |V°| x IVS] matrix, 
where each row corresponds to a vertex of S, and has that vertex’s 
barycentric coordinates on T° as entries. P is defined likewise for 
T!. We then find the function Ô which minimizes Equation 13 as the 
solution to a linear least-squares system, which can be prefactored 
if desired to efficiently transfer many functions. 

We can leverage this formulation to transfer functions from any 
intrinsic triangulation back to the original mesh. In Figure 17, we 
show how this transfer indeed improves the accuracy of PDE solu- 
tions as measured on the original low-quality mesh. This machinery 
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Fig. 18. We construct geodesic paths by flipping edges in a normal co- 
ordinate intrinsic triangulation, as in [Sharp and Crane 2020b]. Normal 
coordinates guarantee a valid path, even on degenerate inputs (left). This 
unlocks advanced applications of geodesics in normal coordinates with the 
same guarantees, such as geodesic loops and Bézier curves (right). 


is enabled because our integer coordinates efficiently and robustly 
compute the common subdivision. More broadly, this paradigm 
opens the door to a wide variety of future finite-element formula- 
tions involving intrinsic triangulations. 


4.5 Flip-Based Geodesic Paths 


The previous sections have demonstrated the value of intrinsic tri- 
angulations as a high-quality basis for discretizing functions on 
surfaces; more broadly, these triangulations also provide simple and 
robust solutions to other tasks across geometry processing. As an 
example, the recent FLIPOUT procedure of Sharp and Crane [2020b] 
computes exact geodesic paths on surfaces via a simple intrinsic 
edge flipping strategy, introducing the geodesic as a path of edges 
in the triangulation. This method is easily implemented in our in- 
teger representation in terms of the mesh operations in Section 3, 
and the resulting geodesic paths may then be recovered with the 
EXTRACTEDGE subroutine. Computing geodesics with our robust 
integer coordinates is particularly appealing, because geodesic algo- 
rithms are notoriously difficult to implement robustly [Sharp and 
Crane 2020b, Section 5.3]. Even the method of Sharp et al. uses the 
signpost data structure, which may fail to reconstruct a connected 
path along the surface for degenerate inputs. In contrast, imple- 
menting FLIPOUT in our integer coordinate representation extends 
the benefits of our approach to this task, including a guarantee 
of valid connectivity in the output (Figure 18, left). It also enables 
higher-level tasks involving geodesic paths to be safely run on low- 
quality input, such constructing geodesic loops on surfaces, and 
even geodesic Bézier curves, using a de Casteljau-style scheme due 
to Morera et al. [2008] as shown in Figure 18, right. 


5 EVALUATION 


We implemented all algorithms in C++; since basic vertex-face adja- 
cency list cannot represent a general A-complex (Section 2.1), we use 
a halfedge data structure for triangle meshes. Timings are measured 
on a single core of an Intel i9-9980XE with 32 GB of RAM. 


Intrinsic Intrinsic 
Delaunay Delaunay 


Method Triangulation Refinement 


Explicit Overlay 100 % - 
Signpost Tracing 96.0 % 69.1 % 
Integer Coordinates 100 % 100 % 


Table 1. The success rate of our method and past approaches for building 
high-quality intrinsic triangulations in the Thingi10k dataset. For each we 
construct a Delaunay triangulation, either on the original vertex set or with 
Delaunay refinement to a 25° minimum angle bound, and attempt to recover 
the connectivity of the common subdivision. The explicit overlay method 
does not support refinement. 


signposts normal coordinates common subdivision 
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Fig. 19. Past methods extracted edges by tracing “signposts” along the mesh, 
which may fail in the presence of degenerate triangles. In contrast, our 
integer coordinates always yield a topologically-valid common subdivision, 
even on extremely poor quality inputs. 


Performance. Generally our data structure is quite fast, computing 
Delaunay refinements for complex meshes in seconds. For example, 
computing the Delaunay refinement Figure 15 takes 0.2s, and the 
Delaunay refinement in Figure 16 (top) takes 0.6s. Because we lazily 
recover intersection geometry from our integer coordinates when 
inserting vertices, routines such as Delaunay refinement which per- 
form many insertions may become moderately expensive on large 
near-degenerate inputs. For instance we take 4 minutes to perform 
Delaunay refinement on 719791 (Figure 20, top) which signposts 
does in 1.5 minutes, but on such meshes signposts generally fails to 
compute a valid common subdivision at the end. Section 6 discusses 
hybrid routines which may give the best of both worlds. 


5.1 Robustness 


We validate robustness by successfully computing Delaunay triangu- 
lations, refinements, and their common subdivisions on all manifold 
meshes in Thingil0k [Zhou and Jacobson 2016]. In particular, we 
used MeshLab to convert each mesh to the PLY file format [Cignoni 
et al. 2008], resulting in 7696 valid manifold meshes. We begin by 
mollifying each mesh to a tolerance of 1077 (Section 3.11). For each 
model we compute the intrinsic Delaunay triangulation (Section 4.1) 
with a tolerance of 1077, as well as an intrinsic Delaunay refinement 
(Section 4.2) with a 25° angle bound. We verify that the algorithms 
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Fig. 20. We fail to compute an explicit mesh of the common subdivision 
following Delaunay refinement on one Thingil0k model (top). Its common 
subdivision would contain 34 million vertices and our program runs out of 
memory. We succeed on a nearly identical model (bottom), whose common 
subdivision contains merely 27 million vertices. 


terminate with the expected conditions. Additionally, we success- 
fully extract an explicit mesh of the common subdivision in both 
cases, except for 1 model in the case of refinement whose common 
subdivision contains around 30 million vertices (Figure 20, top). 
We compare against the explicit overlay representation of Fisher 
et al. [2006] and the signpost representation of Sharp et al. [2019a] 
(Table 1). The overlay representation similarly offers a guarantee of 
valid connectivity, but does not provide a constant-time edge flip 
operation (like normal coordinates do). More importantly it does 
not support operations beyond edge flips and thus cannot perform 
Delaunay refinement. Signposts support a wide range of operations, 
but may not successfully recover the common subdivision on degen- 
erate inputs (Figure 19). The statistic reported here differs from the 
result in Sharp et al. [2019a], because no preprocessing of meshes 
is performed. For refinement Sharp et al. [2019a] do not treat the 
boundary case, so we evaluate only on models without boundary. 


6 LIMITATIONS AND FUTURE WORK 


Limitations. The common subdivisions that we compute after 
Delaunay refinement can be quite large: the mean increase in |V| is 
20x, and the 95" percentile increase is 45x. We emphasize again that 
the common subdivision is not generally a high-quality mesh any- 
way: one should perform numerical computations on the intrinsic 
triangulation instead. The intrinsic mesh has much higher element 
quality and is generally much smaller with a mean increase in |V| of 
3.7x and 95" percentile increase of 7.8x. However, for applications 


that rely on the common subdivision (e.g. Section 4.4), it would still 
be beneficial to explore strategies for simplifying S. 

A related issue is that Delaunay refinement sometimes generates 
meshes with many small triangles. One can prove that Delaunay 
refinement in the plane produces well-graded meshes, meaning es- 
sentially that it only places small triangles in regions with small 
features, and our Delaunay refinement on surfaces seems to behave 
similarly. Nonetheless, on poorly-conditioned input meshes, De- 
launay refinement can insert many small triangles. This can cause 
problems for diffusion-based algorithms (e.g. the logarithmic map 
computation in Figure 16), which use the mean edge length to de- 
termine a suitable diffusion time. We found that computing the 
diffusion time on the original mesh and then performing diffusion 
on the intrinsic triangulation produced the best results. 


Hybrid Data Structures. At this point, there are several intrinsic 
triangulation data structures, but no single one is perfect: 


e Overlay (explicit) provides exact connectivity; flipping can 
be slow; no vertex insertion. 

e Signposts (implicit) provide inexact connectivity; flipping and 
vertex insertion are both fast. 

e Integer coordinates (implicit) provide exact connectivity; flip- 
ping is fast; vertex insertion can be slow. 


We propose a good way to get the best of all worlds would be to 
use a hybrid signpost + integer coordinate data structure. This is 
fully implicit, so you don’t pay the O(n”) cost when you have O(n) 
edges crossing O(n) edges. But, flipping and insertion are both fast, 
and connectivity is exact, if you accept the inserted locations. 

Even further in the implicit direction, storing edge lengths is an 
“optimization” in our data structure. One could just store the nor- 
mal coordinates and original triangulation, recovering edge length 
whenever necessary via a layout operation. This is appealing, since 
it is truly an integer-only representation for intrinsic triangulations. 


General geodesic curves. It would also be natural to use this ma- 
chinery as a representation for general geodesic curves on surfaces, 
which commonly arise in geometry processing tasks such as cutting, 
segmentation, etc.. 
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A PSEUDOCODE 


We assume all algorithms have access to triangulations T° and 
Tt, their edge lengths £9, £1, the normal coordinates n, and the 
roundabouts r. 


Fig. 21. A curve entering triangle jik along edge ij can proceed in 3 ways: 
it can exit along edge ik, in which case it is counted by ae (left); it can exit 
along edge kj, in which case it is counted by cik (center); or it can terminate 
at vertex k(right). This forms the core of algorithm TraceFrom [Reproduced 
from Figure 8 for convenience]. 


Fig. 22. In procedure ExTRACTCURVE, we compute barycentric coordinates 
by laying out a triangle strip in the plane [Reproduced from Figure 9 for 
convenience]. 


Algorithm 1 TraceFrom(¢) 


Input: Any combinatorial crossing ¢ = (ay, p) along some curve y 
lying along TŻ. 
Output: The half of the curve y as a sequence of points 
(G0, 01,.--, Gn, k) along M, where ¢ = ġo and k is the vertex 
at which y terminates. 


1: currentHalfedge — ij 

2: y — [(currentHalfedge, p) ] 

3: while True do 

4: 

5: Tj < currentHalfedge 

6: k — OPPosITEVERTEX(TWIN(currentHalfedge)) 
7z  ifp< ofl then i 

8: currentHalfedge — ik 

9: pep 

10: APPEND(y, (currentHalfedge, p)) 
1: else if p > nij — cik then 

12: currentHalfedge — kj 

13: P— Nt p-nij 

14: APPEND(y, (currentHalfedge, p)) 
15: else 

16: return (y, k) 


Algorithm 2 EXTRACTCURVE(Č) 


Input: Any combinatorial crossing ¢ = (ij, p) along some curve y 
Output: The entire trajectory of y as a sequence of geometric cross- 
ings (a, Zo, Z1,..., Zk, b) along M. 


1: 
2 
3 
4 
5: 
6 
7 
8 
9: 


10: 
11: 
12: 
13: 
14: 


15: 


if runs along edge then 


return e 


: else 


Yfront: b — TRACEFROM() 
Yback, @ — REVERSE(TRACEFROM(C)) 
YCombinatorial S APPEND(Yback: Yfront) 


H — LAYOUTTRIANGLESTRIP(YCombinatorial) 
YGeometric ®© [] 
for ¢ = (ij, p) € YCombinatorial dO 


U, 0 <— INTERSECTIONBARYCENTRIC( Ha, Hp, Hi» Hj) 
z e (ij,p,u,v) 
APPEND(YGeometric z) 


return (a, YGeometric> b), YCombinatorial 


Algorithm 3 ExTRACTEDGE(ab) 


Input: An edge ab € E? 
Output: The entire trajectory of ab as a sequence of geometric cross- 


ings (a, Zo, Z1s+ 2+ Zi b) along M 


1: k — local index of ab about vertex a 


2: ai — argmax7i{rai : 


rai <k} 


3: if rg; = k and nai = —1 then 


4: return ai 

5: else 

6 pock-rg-1 

7: ye ExtractCurvE(NEXT(di), p) 

8 

9 while y does not end at a vertex of V? do 
10: i — endpoint of y 

11: 

12: ¢ — other crossing emanating from j 
13: (i, Z1,...,2Zs, J) — ExTRACTCURVE(¢) 
14: 

15: ADJUSTBARYCENTRICCOORDINATES(Z},..., Zs) 
16: APPEND(Y, (Z1,...,Zs,J)) 

17: return y 


Algorithm 4 SpritFAcE_CAsE1 (ijk, u) 


Input: The location to insert a vertex on T!, as barycentric coordi- 


nates u in a face ijk € F!. 


Output: An updated integer coordinate intrinsic triangulation. 


w y 


: for [= (ij, p) E€ COMBINATORIALCROSSINGS(ijk) do 


zīj [p + 1] — ExrracrGEOMETRICCROSSING(C) 


5: for Jk € CornersOr(ijk) do 


Fig. 23. In SPLITFACE, we do a sequence of line-side tests to compute a value 
of v at each corner. 


T ae p BOY are 


Ci = fe :0<&< Ju € TRIANGLE (i zpi ze leD} 


vi — min (clk UCi) 


jk 
Oj c -Yi 


: if oj = oj, op then 


Oj, Ok — 0,0 


: else if oj > oç, oi then 


Ok, Oi — 0,0 


: else if of = oj,0; then 


i, oj — 0,0 


: for Jk € CornersOr(ijk) do 


Npi — Vit Oj + OK 


: r — UPDATEROUNDABOUTS(Npi, Npj; npk: r) 


ele gl 


pe pp fok — UppaTEEDGELENGTHS(f!, u) 


: R e REGIONFROMNORMALCOORDINATES(Npj, Npj; npk) 
: dp <— RECOVERBARYCENTRIC(R, Z, u) 


Algorithm 5 SpLitFAcE_CAsE2(ijk, u) 


Input: The location to insert a vertex on T!, as barycentric coordi- 


nates u in a face ijk € F!. In Case 1, ijk must be oriented 
so that nij > njk + NK; 


Output: An updated integer coordinate intrinsic triangulation. 


10: 
11: 


for ¢ = (ij, p) € COMBINATORIALCROSSINGS (ijk) do 
zij |p + 1] — ExTRACTGEOMETRICCROSSING(C) 


for JK € CornersOr (ijk) do 


Ci fe :0<&< JK u € Taranote (i, zi [E], zR (2) 
vi — min (c/* UCi) 
jk 


Oj c -Yi 


if o; > o; then 


eS p 


Oj — 0 
: else if oj > oj then 
oj — 0 


: ifv; < clk then 
ij 
k 
: else if vj < oH then 
ij 
k 


vj = vj+e 


Vi — vj+e 
: else 


Eke |g :0<&< eju € TRIANGLE (izl +c) 
€k — min (ed U Ex) 

Vi — Vit EK 

Weer vj rele — ex 

: for Jk € CornersOr(ijk) do 
Npi — Vit Oj + OK 


: r — UPDATEROUNDABOUTS(Npi, Npj; npk» r) 
1 pl pl 
fpi fpj’ fpk 
:ReEe REGIONFROMNORMALCOORDINATES(Npi, Npj: npk) 


: qp «— RECOVERBARYCENTRIC(R, Z, u) 


«— UppaTEEDGELENGTHS(f!, u) 


Al; 


gorithm 6 SPLITEDGE( ij, u) 


Input: The location to insert a vertex on Tt, as barycentric coordi- 


nates u on a halfedge ij c H! 


Output: An updated integer coordinate intrinsic triangulation 


1 


: if njj 2 0 then 
k «— OPPOSITEVERTEX( ij) 


3: SpLITFACE(ijk, (ui, uj, 0)) 
4: FLIPEDGE(i/) 
5: else 
6: ke OpposITEVERTEX( ij) 
7: if ININTERIOR(i/) then 
8: | — OpposITEVERTEX( ji) 
9: T! e insert a vertex p along ij 
10: pj, pi — Nij, Mij 
11: npk — Max(Ngj, NjK, 0) 
12: 
13: r — UPDATEROUNDABOUTS (pi, Np j; npk: Nph r) 
14: ti Op oe bi «— UPDATEEDGELENGTHS (£1, u) 
15: qp = uig? + ujjai 
Algorithm 7 REMOVEVERTEX(i) 


Input: An inserted vertex i 


Output: Updated triangulation i removed 


1: 
2: 
3: 


4 


while i has degree > 3 do 
ij < flippable edge incident on i 
FLIPEDGE(i/) 


: DELETEVERTEXANDINCIDENTEDGES (i) 


Algorithm 8 CompuTECOMMONSUBDIVISION() 


Input: Nothing beyond the usual data (i.e. T°, T,...) 
: T? « index common subdivision vertices 


2 

3: 

4: polygons + [] 
5: for ijk € F! do 
6 if e) = 0 then 

7 APPEND(polygons, SUBDIVIDEFACE_CAasE1(ijk, Z”)) 
8: else 

9 APPEND(polygons, SUBDIVIDEFACE_CasE2(ijk, Z”)) 


1: for i € V! do 


2: ke H 

33 Qe — q 

4: Q} É i (i, 1) 

5: for ab € E° do 

6 y = (a, Z1, ..., Z}, b) — EXTRACTEDGE(ab) 
7: for z = (ij, p,u,v) € y do 

8 ke Ti lp +1] 

9 Q? — (ab, u) 

20: Q} — (ij,v) 


21: return polygons, Q°, Q! 


Algorithm 9 DELAUNAYREFINEMENT(@min ) 


Input: A minimum allowed angle Onin. 


Output: An intrinsic triangulation T! whose corner angles are all at 


least Omin 
1: FrrPTODELAUNAY() 
2: while T! has triangles with angles less than Omin do 
3: ijk — any triangle with an angle less than Onin 


Ue < circumcenter barycentric coordinates 

uc — Ue — (1/3, 1/3, 1/3) 

V — BARYCENTRICOFFSETTOTANGENTVECTOR( due) 

c — Exp(BARYCENTER(ijk), V) 

if c lies inside the mesh then 
INSERTCIRCUMCENTER(ijk) 

else 


lm <— boundary edge separating c from ijk 
p <— SpLiTEDGE(Im, 0.5) 


FLrpTODELAUNAY() 


PO T ie SY cae, OR SRS sh eS NO) Ger a GS OB a 


Y 
2 


ball = {i € V! : DyKsTRADISTANCE(E!, i, p) < lim} 
for i € ball do 
REMOVEVERTEX(i) 


N N 
N e 


23: FrrpPTODELAUNAY() 


B BARYCENTRIC COORDINATES RECOVERY 


In Section 3.4, we recover the barycentric coordinates of a newly 
inserted vertex on T° via interpolation along a polygonal subre- 
gion R of triangle abc € F°. Here we give a full expression for the 
necessary small linear system. 

Precisely, let 3 < p < 6 denote the number of corners of R. Let the 


m™ corner of R have barycentric coordinates u™ „u m) y 


abc € F’ and barycentric coordinates o, o, o™ 
all of which are know. We also know the barycentric coordinates v; 
for p in ijk. We then want to solve for the corresponding ua on abc. 
We proceed in two steps: first, we express v as a linear combination 
E of the v0”), Then, we apply this same linear combination to the 
u\™) to obtain u. Concretely, we first solve for the minimum-norm 


solution of the underdetermined system 


ROMER) RONEL 
i i & Ui 
(0) yi) vP = llo; (14) 
lo to tal: a 
t U 
ra) Uk Uk 7 k 


and then set 


Ua := >» ui) Em: Up := > age Ems 
m 


w > us Gas (15) 


Note that while one often seeks a nonnegative č, any solution will 
suffice here: we only use ë to interpolate in Equation 15. 


C DELAUNAY REFINEMENT DETAILS 
C.1 Removing Extra Vertices 


When Chew’s second algorithm splits an edge, it removes all inserted 
circumcenters within a geodesic ball centered at the edge’s midpoint. 
These vertices must be removed, but it is okay to removes additional 
interior inserted vertices. Shewchuk [1997, Section 3.4.2] observes 
that the algorithm can only perform finitely many edge splits. As 
long as one removes all interior inserted vertices within the geodesic 
ball—and never removes vertices along the boundary—the algorithm 
will still perform only finitely many edge splits. Hence, it must 
terminate as usual following the final edge split, even if one removes 
extra circumcenters during edge splits. 


C.2 Proof of Correctness on Watertight Meshes 


Here we seek to prove that DELAUNAYREFINEMENT (Algorithm 9) 
succeeds, in the basic case of a closed surface with bounded cone 
angles. We will not prove the more general boundary case here, but 
experimentally we observe success on a large dataset (Section 5.1). 


THEOREM C.1 (DELAUNAY REFINEMENT, NO BOUNDARY). On meshes 
without boundary, with vertex angle sums at least 60°, Algorithm 9 
produces a Delaunay mesh with triangle corner angles at least 30°. 


Proof. By definition, DELAUNAYREFINEMENT only terminates 
when the triangulation is a Delaunay triangulation which satisfies 
the angle bound, so we just need to prove that termination occurs 
after a finite number of iterations. We will show this by establishing 
that DELAUNAYREFINEMENT maintains a minimum spacing between 


all vertices in the mesh, so the number of insertions is bounded by 
surface area. Our argument will generally follow the planar proof 
of Shewchuk [1997, Section 3.2.1], though extra care is needed in 
the intrinsic case, where self edges may connect a vertex to itself. 

In particular, we consider the length of the shortest edge in the 
initial mesh’s intrinsic Delaunay triangulation, ô := minj;; fij. We 
will show that the minimum edge length in each subsequent Delau- 
nay triangulations is at least ô. Then all vertices must be separated 
by a distance at least 6, since Lemma C.3, each vertex is connected 
to its geodesic nearest neighbor. Hence, each vertex is contained in 
an open disk of radius 55 which is disjoint from all other disks. As 
the input mesh has finite surface area, we conclude that Algorithm 9 
can only insert finitely many vertices, and thus must terminate. 

It remains to show that DELAUNAYREFINEMENT never creates an 
edge of length less than 6. It is convenient to convert the angle 
bound @ to a circumradius-to-shortest-edge ratio bound B = se 
[Shewchuk 1997, Section 3.1]. Having corner angles at least a = 30°, 
is equivalent to a circumradius-to-shortest-edge ratio of at most 
B = 1, and thus we insert the circumcenters of triangles with B > 1. 

We proceed by induction. All initial edges have length at least 6 
by definition. Now consider inserting vertex i at the circumcenter of 
triangle jkl with circumradius R. Since we only split triangles with 
B > 1, and jkl’s edges have length at least ô, we must have R > ô. 
By Lemma C.4 all new edges in the Delaunay triangulation must 
be incident on i, and since jkl had an empty geodesic circumcircle, 
there can be no other vertices within distance R > 6. Thus new all 
edges to other vertices have length at least 6. We must now consider 
self edges connecting the new vertex i to itself. 

Gluing together the two ends of a self edge yields a loop; we will 
split into cases based on the homotopy class of this loop on the 
punctured surface (before the insertion of i). First, note that the loop 
cannot be contractible to a point, since the original edge is geodesic. 
Then we will split in to two cases: either the loop contracts around 
a single vertex, or it does not. 

If the loop contracts around a single vertex, then 
the self edge encloses a degree-1 vertex. The degree- 
1 vertex must have distance at least R to the inserted 
vertex, and has angle sum at least 60°. Thus, by the 
law of cosines, the length of the self edge must be 
at least 


VR2 + R2 — 2R? cos 0 = RY2(1 — cos 0). 


Since cos 60° = 3 and 1 — cos @ is increasing with 
0, this shows that the self edge has length at least 
R whenever @ is at least 60°. 

If the loop is not in a homotopy class contractible 
about a single vertex, then the shortest loop ymin 
in the homotopy class is non-constant. By Lemma C.5, we can take 
Ymin to touch some vertex a, and note that since ymin is the shortest 
loop our original self edge must be at least as long as ymin. Then by 
Lemma C.3 a has an edge at least a long as ymin, and thus the self 
edge has length > |ymin| 2 ô. 

Thus, we conclude that Algorithm 9 never introduces an edge 
of length less than ô, which means that it must terminate after 
inserting finitely many vertices. o 


>42 a —cos 0) 


Lemma C.2. For any pair of vertices i, j € V, let Tj; be the set of 
non-constant geodesics connecting i to j. Then 


dij = inf length(y) > 0. 
yeliy 


Proor. This follows directly from [Indermitte et al. 2001, Proposi- 
tion 1], which states that for any L > 0, the number of geodesic arcs 
from i to j of length at most L is finite. Since any geodesic of length 
0 is constant, and thus not in Ij;, this implies that dj; > 0. o 


LEMMA C.3. For any vertex i € V, the intrinsic Delaunay triangu- 
lation contains an edge to i’s nearest neighbor. 


Proor. This is a standard result, which we include for complete- 
ness. Let j be i’s nearest neighbor, i.e. j := argmin ‚dij. Note that 
j may equal i, and di; > 0 by Lemma C.2. Consider the disk D of 
radius dj; centered at i. Since j is i’s nearest neighbor, D contains no 
vertex other than i. Thus, the circle which goes through i and j and 
is tangent to D at j has empty interior, and its boundary contains no 
vertices other than i and j. We conclude that ij is in the Delaunay 
triangulation [Bobenko and Springborn 2007, Definition 3]. oO 


Lemma C.4. All edges created in DELAUNAYREFINEMENT following 
the insertion of a vertex i and flipping to Delaunay are incident on i. 


Proor. Again, we follow the planar proof of Shewchuk [1997, 
Lemma 12]. We wish to prove that all Delaunay edges which are not 
incident on i were Delaunay before inserting i. This follows from 
the fact that edges of a Delaunay triangulation satisfy an empty 
circumcircle condition [Bobenko and Springborn 2007, Definition 
3]. If an edge’s circumcircle is empty after inserting vertex i, it must 


have been empty before too, so the edge was already Delaunay. o 


Lemma C.5. Any geodesic loop y is isotopic to a geodesic loop y’ of 
the same length which touches a vertex. 


Proor. y can “slide” until it touches a vertex without changing 
its length. Precisely, consider a unit-speed motion of y within the 
surface along its outward normal direction. During the motion, 
4 lyl= f k(s) ds, where x is the geodesic curvature of y. Since y is 
a geodesic, x = 0: its length does not change. Thus we can construct 


y’ by sliding y along the surface until it touches a vertex. o 


As an aside, we note that geodesic loops which do not touch a 
vertex only occur in non-generic configurations. 


Je 


Fig. 24. The relevant angles for Theorem D.1. 


D  SIMPLICIAL VERTEX REMOVAL 


In Section 3.6, we consider removing a vertex by flipping edges until 
the vertex has degree three and then deleting it. Past work has also 
proposed this approach in the purely-topological setting [Schaefer 
et al. 2002, Section 5.4], but here we must respect geometric con- 
straints. In particular, edges can only be flipped geometrically if they 
are contained in a convex quadrilateral (Section 3.3). Here we prove 
that flipping edges to remove vertices is indeed a viable strategy in 
the Euclidean setting as well: one can always find an edge to flip. 


THEOREM D.1 (VERTEX REMOVAL, SIMPLICIAL). If a vertex i in a 
simplicial complex has cone angle 2x and degree d > 3, then some 
edge ij incident on i can be flipped to decrease the degree of i. 


Proor. Recall that an edge can be flipped if both endpoints will 
have degree at least 1 after the flip, and the edge is contained in a 
convex quadrilateral (Section 3.3). As always, the convex quadri- 
lateral is defined in the sense of the intrinsic geometry determined 
by edge lengths. The endpoint degree constraint is automatically 
satisfied on a simplicial complex, so we only need to show that the 
geometric convexity constraint is satisfied, which is equivalent to 
showing that all angles of the edge’s quadrilateral are at most z. 

Denote the neighboring vertices of i as jọ, with j,4, etc. implicitly 
indexed modulo the vertex degree d (Figure 24). The outer angles 
Lijpejg and 2; jpk are corners of Euclidean triangles, and thus are 
necessarily at most 7, so we need to find an edge ij, for which the 
angles Zj,_sijgy, and Zjx,, jy jg are also at most z. 

First we consider the inner corners Z/j,_,ij,,,. At most two of 
these angles can be greater than z. To see why, suppose there were 
three 4j,_,ij,,; > T. Since the degree of i is d > 3, then some pair 
of those three large angles would correspond to disjoint angular 
sectors around the vertex, and summing their angles yields a value 
greater than 22, which is impossible because the angle sum of i is 
2. Thus all but at most two of the edges incident on i have inner 
corners with angle at most z. 

Likewise, at least three of the outer corners Z j,,, j, j,_; are at most 
x. This is because the sum of all d outer corners must be (d — 2)z. 
Since they are nonnegative, at most d — 3 of them can be strictly 
greater than x, implying that at least 3 will be z. 

Thus at least three outer corners are at most z, and at most two 
of the inner corners are not at most x, so there must be at least one 
edge for which both the inner and outer corners are at most 7. This 
edge can then be flipped, reducing the vertex degree. o 


Importantly, this proof does not handle the full general case of a 
A-complex, where there may exist self-edges which cause flips to 
not make progress. However, we note that Sharp and Crane [2020b, 
Appendix A] proves that a similar flip-removal strategy works in the 
case of a A-complex, and we conjecture that an analogous technique 
could be applied to generalize Theorem D.1. Also, note that the 
“equality” case of Theorem D.1 is a possibility, such as a degree 
four cross configuration where all angles = 2/2. Fortunately the 
resulting skinny triangle after the edge is a non-issue, because the 
center vertex is about to be removed. 


